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ABSTRACT 

We use the NIHAO (Numerical Investigation of Hundred Astrophysical Objects) cos¬ 
mological simulations to study the effects of galaxy formation on key properties of 
dark matter (DM) haloes. NIHAO consists of w 90 high-resolution SPH simulations 
that include (metal-line) cooling, star formation, and feedback from massive stars and 
SuperNovae, and cover a wide stellar and halo mass range: 10® < AA/Mq < 10^^ ( 
10®-® < Mhaio/M 0 < 10^®-®). When compared to DM-only simulations, the NIHAO 
haloes have similar shapes at the virial radius, i?vir, but are substantially rounder 
inside ~ O.lAvir- In NIHAO simulations c/a increases with halo mass and integrated 
star formation efficiency, reaching ^ 0.8 at the Milky Way mass (compared to 0.5 in 
DM-only), providing a plausible solution to the long-standing conflict between obser¬ 
vations and DM-only simulations. The radial profile of the phase-space Q parameter 
(p/cr®) is best ht with a single power law in DM-only simulations, but shows a flat¬ 
tening within 0.1i?vir for NIHAO for total masses M > IO^^Mq. Finally, the global 
velocity distribution of DM is similar in both DM-only and NIHAO simulations, but in 
the solar neighborhood, NIHAO galaxies deviate substantially from Maxwellian. The 
distribution is more symmetric, roughly Gaussian, with a peak that shifts to higher 
velocities for Milky Way mass haloes. We provide the distribution parameters which 
can be used for predictions for direct DM detection experiments. Our results underline 
the ability of the galaxy formation processes to modify the properties of dark matter 
haloes. 

Key words: Galaxy: disc, evolution, structure - galaxies: disc, evolution, interactions, 
structure - methods: numerical, N-body simulation 


1 INTRODUCTION 

N-body numerical cosmological simulations have proven 
to be powerful tools for modeling the properties of the 
3-dimensional dark matter distribution in galactic haloes 
(Jing & Suto 2002; Allgood et al. 2006; Bett et al. 2007; 
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Hayashi et al. 2007; Maccio et al. 2007; Neto et al. 2007). 
However, these dark-matter-only cosmological simulations 
only take into account the effects of gravity. Globally, the 
baryons only make up Db/Dm ~ 15% of the mass budget 
(Planck Collaboration et al. 2014), and thus on large scales 
their influence can be reasonably ignored. However, baryons 
can dissipate energy, leading to large mass fractions near 
the centers of galaxies. Additional processes such as the vio- 
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lent explosion of stars as supernovae can have non-negligible 
impact on the dark matter distribution. 

Using the ACDM cosmology, dark-matter-only galaxy 
simulations can make precise predictions of the shape and 
the internal matter and velocity distribution of dark matter 
haloes (Maccio et al. 2008). On average, such dissipation¬ 
less cosmological simulations predict the ratio of the minor 
principal axis to the major to be (c/a) ~ 0.6. Observations, 
however, indicate that haloes are more spherical. Ibata et al. 
(2001) used the Sagittarius Stream (and a Dehnen & Binney 
(1998) halo mass distribntion) to pnt a lower bound on the 
shape of the inner Milky Way halo of c/a 0.8, where the 
halo was defined to be between 20 kpc < r < 60 kpc. Other 
recent work fonnd similarly spherical shapes for the inner 
halo of the Milky Way (Helmi 2004; Martmez-Delgado et al. 
2004). 

In a recent, detailed study. Law & Majewski (2010) pre¬ 
sented a new N-body model (with a logarithm dark mat¬ 
ter halo)^ for the tidal disruption of the Sagittarius galaxy 
which led them to conclude that our Galaxy has an oblate 
potential with a minor to major axis ratio of (c/a)^ = 0.72. 
This latter value can be seen as a upper limit on the shape 
of the matter distribution, which again suggests a possible 
tension between results of pure CDM simulations and local 
measurements of halo shape. 

One possible solution to this problem could be re¬ 
lated to the inclusion of a dissipative component such as 
baryons in cosmological simulations. Katz & Gunn (1991) 
and Dubinski (1994) were among the first to show that 
haloes in dissipative simulations were systematically more 
spherical than corresponding haloes in dissipationless sim¬ 
ulations. Kazantzidis et al. (2004) analyzed the effect of 
baryons on the shape of dark matter haloes using high res¬ 
olution hydrodynamical simulations of clusters of galaxies 
and, again, found a reduced triaxiality in dissipational simu¬ 
lations. That result was confirmed in several further studies 
(Springel et al. 2004; Debattista et al. 2008; Pedrosa et al. 
2009; Tissera et al. 2010). While there is an agreement that 
baryons tend to make dark matter haloes rounder, a quanti¬ 
tative prediction of this effect strongly depends on the pre¬ 
scription used to model baryons in cosmological simulations, 
which is not a trivial task. 

A large fraction of simulated galaxies used for study¬ 
ing the impact of baryons on the DM distribution were ei¬ 
ther plagued by the so called over-cooling problem (which 
produced galaxies that were unrealistically dominated by 
baryons in their central region Kazantzidis et al. 2004) or 
were missing important aspects of the galaxy formation pro¬ 
cess such as star formation and feedback (Abadi et al. 2010). 

A noticeable exception has been the recent work by 
Bryan et al. (2013), which employed the OWLs simulations 
(Schaye et al. 2010), a large scale simulation that produced 
realistic large galaxies and cosmic structure in a full cos¬ 
mological framework. These simulations found that baryons 
substantially changed c/a of massive galaxies. On the other 
hand they were able to resolve galaxies on the scale of the 
Milky Way with only few thousand particles. 

Shape is not the only parameter of dark mat- 


^ An NFW profile will give similar results (Law &; Majewski 
2010) 


ter haloes that baryons might modify. A lot of re¬ 
cent attention has focused on the modification (expan¬ 
sion and core creation) of density profiles (Gnedin et al. 
2004; Abadi et al. 2010; Governato et al. 2010; Maccio et al. 
2012; Pontzen & Governato 2012; Di Gintio et al. 2014b; 
Dutton et al. 2015; Ofiorbe et al. 2015), which helps re¬ 
duce the tension on small-scales between observations (e.g. 
Oh et al. 2015, and references therein) and DM-only pre¬ 
dictions (Toilet et al. 2016; Dutton et al. 2016b). Another 
interesting quantity that is a slight variant on the density 
profile is the so called coarse-grained phase-space density 
Q = p/a^ that takes into account particle velocities as well 
as matter density. DM-only simulations show that radial 
profiles of Q are well fit with a simple power law p/cr® oc r“ 
(Taylor & Navarro 2001). 

Analytical work found a characteristic value for the 
power law slope of the Q profile, a = 1.944, for isotropic 
structures (Austin et al. 2005). This value also reproduces 
realistic radial velocity dispersions through cosmological 
simulations (Dehnen & McLaughlin (2005), but see also 
Schmidt et al. (2008) for a more critical interpretation). It is 
particularly interesting to check if baryons break this sim¬ 
ple power law behavior of the dark matter coarse-grained 
phase-space density in simulated galaxies. 

The dark matter velocity distribution is another very 
important parameter, since it has profound implications for 
the predicted dark matter - nucleon scattering rates in di¬ 
rect detection experiments (Kuhlen et al. 2010). Most mod¬ 
els assume a Maxwell-Boltzmann (MB) velocity distribution 
function. Therefore, a departure from the MB distribution 
might change the predicted rate of events for dark matter 
models in which the scattering is sensitive to the high veloc¬ 
ity tail of the distribution, and for experiments that require 
high energy recoil events (e.g. Vergados et al. 2008). 

Recent high resolution DM-only simulations have re¬ 
ported substantial departures from a Maxwellian shape 
(Vogelsberger et al. 2009; Kuhlen et al. 2010), with a deficit 
near the peak and excess particles at high speeds. In con¬ 
trast, an analysis of the ERIS hydrodynamical simulation of 
a single Milky Way like galaxy, showed a different behav¬ 
ior. Its velocity function maintains the Maxwellian distribu¬ 
tion and shows a greater deficit than MB at high velocities 
(Pillepich et al. 2014). 

In this paper we revisit the issues of the effect of baryons 
on dark matter properties using the galaxies from the NI- 
HAO project (Wang et al. 2015). The NIHAO project is a 
large suite of high resolution simulated galaxies. The sam¬ 
ple includes more than 90 high resolution zoom regions with 
halo masses ranging from ~ 10^° to IO^^Mq. Each high 
resolution galaxy is resolved with at least 600, 000 elements 
(DM-I-GAS) and up to several millions [see Wang et al. 
(2015) for details]. It is the largest sample of high resolution 
galaxies to date. The simulations show remarkable agree¬ 
ment with the stellar mass-halo mass relationship across 
five orders of magnitude of stellar mass. Thus a fully sam¬ 
pled volume would follow the observed galaxy stellar mass 
function. The simulations thus offer a unique tool to study 
the modifications that baryons induce to dark matter since 
it simultaneously combines high spatial and mass resolution 
with a statistical sample of galaxies across three orders of 
magnitude in halo mass. 

We have decided to focus this paper on three main prop- 
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erties: the halo shape, since it is a widely studied halo prop¬ 
erty and it is possible to directly measure it, especially in 
our own galaxy. The (pseudo) phase-space density, since it 
has been claimed to be an universal property of DM haloes 
(Taylor & Navarro 2001) and it is somehow complementary 
to the DM density profiles which have been discussed in the 
NIHAO IV paper, and the DM velocity distribution that has 
important consequences on the predicted rate of events in di¬ 
rect detection experiments. Another interesting property is 
the angular momentum distribution of the DM, but this will 
be the topic of a different NIHAO paper more closely linked 
to the process of disc formation (Obreja et al. in prep.) 

This paper is organized as follows: in Section 2 we dis¬ 
cuss the simulations including the feedback prescription used 
in creating the NIHAO galaxies; In Section 3 we discuss in 
detail the effects of baryons on the dark matter halo shape, 
pseudo phase space density and the velocity distribution. 
Section 4 offers a summary of our results, as well as their 
impact. 


2 SIMULATIONS 

We study simulations from the NIHAO project (Wang et al. 
2015), that use baryonic physics from an updated ver¬ 
sion of the MaGICC simulations (Stinson et al. 2013). The 
initial conditions are created using the latest compila¬ 
tion of cosmological parameters from the Planck satellite 
(Planck Collaboration et al. 2014); namely Qa ~ 0.6825, 
Dm =0.3175, Ho = 67.1, ug = 0.8344, n = 0.9624 and 
Db = 0.0490. The haloes to be re-simulated with baryons 
and higher resolution were extracted from three different 
pure N-body simulations with box sizes of 60, 20 and 15 
/i“^Mpc, more information on the DM-only simulations can 
be found in Dutton & Maccio (2014). 

The NIHAO project was designed to study galaxy for¬ 
mation over a wide mass range, from dwarf galaxies to mas¬ 
sive spirals like the Milky-Way. The simulations maintain 
the same numerical resolution across the entire mass range. 
This means that all haloes are resolved with at least 600,000 
elements (DM-l-stars+GAS) within the virial radius, up to 
several millions. This fixes numerical resolution has been 
achieved by varying both the size of the initial box from 
which the haloes to be zoomed have been chosen, and the 
number of zoom levels of each simulation (see figure 2 in 
Wang et al. (2015)). 

The zoomed initial conditions were created using a mod¬ 
ified version of GRAFIC2 (Bertschinger 2001; Penzo et al. 
2014). All haloes/galaxies presented in this paper are the 
most massive object in their respective high resolution vol¬ 
ume, in other words we will only present results, at all 
masses, for central haloes. The starting redshift is Zstart = 
100, and each halo is initially simulated at high resolution 
with DM-only using pkdgrav (Stadel 2001). More details 
on the sample selection can be found in Wang et al. (2015). 

We refer to simulations with baryons as the Hydro sim¬ 
ulations or NIHAO, while we will use the term N-body or 
DM-only for the collisionless simulations. 

The hydro simulations are evolved using an improved 
version of the SPH code gasoline (Wadsley et al. 2004). 
The code includes a subgrid model for turbulent mixing of 
metals and energy (Wadsley et al. 2008), heating and cool¬ 



FigUre 1. Stellar mass - halo mass relation for the NIHAO galaxies used 
in this work. All simulations have more than 600,000 particles in their virial 
radius (see Wang et al. 2015 for more details). The solid lines represent the 
most commonly used abundance matching results (see text). 

ing include photoelectric heating of dust grains, ultraviolet 
(UV) heating and ionization and cooling due to hydrogen, 
helium and metals (Shen et al. 2010). 

For the NIHAO simulations we have used a revised 
treatment of hydrodynamics described in Keller et al. (2014) 
that we refer to as esf-Gasoline2. Most important is the 
Ritchie & Thomas (2001) force expression that improves 
mixing and shortens the destruction time for cold blobs (see 
Agertz et al. (2007)). esf-Gasoline2 also includes the time- 
step limiter suggested by Saitoh & Makino (2009), which is 
important in the presence of strong shocks and temperature 
jumps. We also increased the number of neighbor particles 
used in the calculation of the smoothed hydrodynamic prop¬ 
erties from 32 to 50. 

2.1 Star formation and feedback 

The star formation and feedback modeling follows what was 
used in the MaGICG simulations (Stinson et al. 2013). Gas 
can form stars when it satisfies a temperature and a den¬ 
sity threshold: T < 15000 K and nth > 10.3 cm“®. Stars 
can feed energy back into the ISM via blast-wave super¬ 
nova (SN) feedback (Stinson et al. 2006) and via ionizing 
radiation from massive stars before they turn in SN. Met¬ 
als are produced by type H SN, type la SN. These, along 
with stellar winds from asymptotic giant branch stars also 
return mass to the ISM. The metals affect the cooling func¬ 
tion (Shen et al. 2010) and diffuse between gas particles 
(Wadsley et al. 2008). The fraction of stellar mass that re¬ 
sults in SN and winds is determined using the Chabrier 
(2003) stellar Initial Mass Function (IMF). 

There are two small changes from the MaGIGG sim¬ 
ulations. The change in number of neighbors and the new 
combination of softening length and particle mass means 
the threshold for star formation increased from 9.3 to 10.3 
cm“^, The increased hydrodynamic mixing necessitated a 
small increase of pre-SN feedback efficiency, cesf, from 0.1 
to 0.13. This energy is ejected as thermal energy into the 
surrounding gas, which does not have its cooling disabled. 
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Table 1. Properties of the four selected “test” galaxies. M* is the total stellar mass within .Rvir- Note that in the NIHAO project the 
“name” of a galaxy indicates the mass of the halo in the DM-only simulation. 


Simulation 

name 

M* 

[Mo] 

-^vir 

[kpc] 

Softening 

[pc] 

rriDM 

[Mq] 

TTlgas 

[Mo] 

]\fDM 

vir 

N* 

vir 

Mass bin 

name 

Mass range 

[Mo] 

gl.92el2 

1.59e+ll 

281 

931 

1.74e+06 

3.16e+05 

1,200,667 

2,467,821 

M3 

7.5e-|-ll< Mvu <3.5e+12 

g5.02ell 

1.46e+10 

179 

465 

2.17e+05 

3.95e+04 

2,408,247 

1,821,403 

M2 

2.5e-|-ll< Mvir <7.5e+ll 

gl.OSell 

8.47e+08 

105 

465 

2.17e+05 

3.95e+04 

520,108 

108,025 

Ml 

7.5e-|-10< Mvir <2.5e+ll 

g4.99el0 

1.24e+08 

78 

207 

1.90e+04 

3.48e+03 

732,808 

52,511 

MO 

4.0e-|-09< Mvir <7.5e+10 



Figure 2. Histogram of the ratio between stellar half mass radius, Ri/ 2 ^ 
and the viral radius, itvir- For reference, the three radii at which we will 
compute the halo shape are: -Rvir; 12% -Ryir (red line) and 5% -Ryir (blue 
line). The leftmost histogram shows the distribution of the stellar particle 
softening through NIHAO. 

Most of this energy is instantaneously radiated away, and 
the effective coupling is of the order of 1%. 

2.2 The galaxy sample 

For this project we use a total of 93 galaxies from the NI¬ 
HAO sample, all galaxies are “centrals” in their respective 
haloes and are resolved with have more than 600,000 parti¬ 
cles (dm+gas-|-stars) within the virial radius.^ Since our aim 
is to study the impact of galaxy formation on the dark mat¬ 
ter distribution it is very important to use realistic simulated 
galaxies, i.e. galaxies able to reproduce the observed scaling 
relations. One of the key success of the NIHAO galaxies is 
to be able to reproduce the observed stellar mass-halo mass 
relation across the whole mass spectrum, which extends for 
more than 5 orders of magnitude in stellar mass, as shown in 
Fig. 1. As detailed in Wang et al. (2015) the NIHAO galax¬ 
ies are also able to reproduce the the stellar mass - halo mass 
relation also at higher redshift, and have realistic star forma¬ 
tion rates for their stellar masses. NIHAO galaxies are also 
consistent with the observed gas content of galaxy discs and 

^ In this paper we have used the latest compilation of NIHAO 
galaxies, including several new “central” galaxies that were com¬ 
pleted after the publication of the first NIHAO paper. 


haloes (Stinson et al. 2015; Wang et al. 2016; Gutcke et al. 
2016). Thanks to the unprecedented combination of high res¬ 
olution and large statistical sample, the NIHAO suite offers 
a unique tool to study the distribution response of several 
DM properties to galaxy formation. 

Since among our goals there is the study of the shape 
of the DM distribution at different radii, it is interesting 
to check how these radii compare to the size of the central 
stellar region. Fig. 2 shows a histogram of the ratio between 
the stellar half-mass radius, R 1 / 2 , and the viral radius, i?vir. 
The three key radii at which we will compute the halo shape 
are: i?vir, 12% i?vir and 5% Rvii- As the plot shows all these 
radii are larger than the size of the stellar body and are well 
above the softening of the simulations which is shown by the 
leftmost histogram in the figure (see also Table 1 and Wang 
et al. 2015). Finally we have decided to use as minimum 
radius for our profiles 2% of i?vir, which again is well within 
the softening of all simulations. A more detailed discussion 
of the “convergence radius” of our simulations can be found 
in Toilet et al. (2016). 

In order to organize our results we have decided to show 
detailed results (as for example the triaxiality parameter as a 
function of radius) for only four “test” galaxies. These galax¬ 
ies have masses equally spaced by half a dex ranging from 
5 X 10^° to 10^^ and are listed in Table 1. A visual impression 
of our galaxies is shown in figure 3, and was created using the 
Monte Carlo radiative transfer code SUNRISE (Jonsson 2006). 
The image brightness and contrast are scaled using arcsinh 
as described in Lupton et al. (2004). Around the mass each 
of those galaxies we have then constructed a (total) mass 
bin of size of 0.5 dex (i.e. the mass bin around gl.OSell goes 
from 7.5 x 10^° to 2.5 x IO^^Mq). No other selection criterion 
is used besides the mass. In the following, for all inspected 
DM properties (shape, phase-space and velocity distribu¬ 
tion) we will present results for our single test galaxies and 
then average results for the 4 mass bins: M0-M3, listed in 
increasing mass order. 


3 RESULTS 

To determine the properties of the dark matter halo we con¬ 
sidered all dark matter particles within Rvir (the radius that 
encompasses a density equal to 200 times the critical density 
of the Universe) found using the Amiga Halo Finder (AHF) 
(Knollmann & Knebe 2011). The study makes a direct com¬ 
parison between the haloes that form in the DM-only simu¬ 
lations (black symbols in each Fig.) and those that form in 
the hydro simulations (red symbols). 

The comparison is made between the axis ratios {b/a, 
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Figlir© 3. Face-on (upper panels) and edge-on (lower panels) views of our four selected “test” galaxies: g4.99el0, gl.OSell, g5.02ell, gl.92el2 from left to 
right. Galaxies have been processed through the Monte Carlo radiative transfer code SUNRISE. Images are 50 kpc on a side. 


c/a and T the triaxiality), the pseudo phase space density 
(Q), and the velocity distribution. We investigate the veloc¬ 
ity dispersion both globally and in the Solar neighborhood. 

3.1 Halo shape 

For each simulation (DM-only and hydro), we define the 
principal axes of the shape ellipsoid such that a > b > c and 
define the ratio parameters, as s = c/a, q = b/a,p = c/b. To 
calculate the shape of the dark matter within a given radius, 
we begin with the assumption of a spherical ellipsoid, such 
that a = b = c. We then compute the inertia tensor pj, 
defined as: 

lij — Xj /, 

where m is the particle mass, a is the particle index and i,j 
refer to the coordinates. The radius is defined to be = 
Xa + vh/ff + Da / (Kazantzidis et al. 2004). 

The eigenvalues of this matrix produce new values for 
s,q,p. We iterate, using the eigenvalues of the previous ma¬ 
trix as the new assumptions of s and q, until the fractional 
difference of s,q,p converges to a tolerance value of 10~® 
(Maccio et al. 2008). 

Finally we define the triaxiality parameter as T = [1 — 
(&/a)^]/[l — (c/a)^]. A prolate halo has T — 1, and oblate 
halo has T = 0, while a triaxial halo has T ^ 0.5. 

In Fig. 4, we present the relation between the minor to 
major axis ratio (c/a, left) the middle to major [b/a, center) 
and the triaxiality parameter (right) as a function of radius. 
The solid line is for the integral measurement of the shape 
(i.e. c/a{< r)) while the dotted line is for the differential 
one. In each panel the last point marks the virial radius of 
the halo. As discussed previously for clarity and brevity we 
present results only for our 4 test galaxies (see Table 1) out 
of 93 in our sample. 


The trend is similar across all masses: baryonic physics 
processes such as cooling and stellar feedback modify the 
shape of the DM distribution. The effect of baryons physics 
is strongest at the center of the halo. In general, the hydro 
simulations have a less varying shape as a function of radius 
w.r.t Nboby sims and the hydro shape is always rounder 
than the dark matter. At small radii, the dark matter be¬ 
comes more triaxial. Near i?vir, the shapes of the hydro and 
DM-only simulations converge to their most spherical shape 
in all the galaxies except the most massive one, gl.92el2. 

Fig. 5 shows the same quantities as in Fig. 4 but this 
time averaged over the four mass bins M0-M3 (see Table 1), 
with the grey area representing the la scatter around the 
mean. The radius at which DM-only and Hydro results de¬ 
part from each other moves to larger radii as the mass grows. 
In the lowest mass bin (MO, lowest row) there an apprecia¬ 
ble difference between DM-only and Hydro only inside few 
kpc, while in the most massive bin (M3, uppermost row) the 
difference is already substantial at 100 kpc. This is true for 
all three diagnostics. Overall the general trends are quite 
in agreement with the ones from our test galaxies, which 
appear to provide a good representation of the general pop¬ 
ulation, on the other hand this figure also points out the not 
so negligible galaxy-to-galaxy variations in the response of 
haloes to baryons. 

Debattista et al. (2008) described how the condensation 
of baryons into galaxy centers affect the dark matter shape. 
DM-only simulations produce haloes with a prolate shape, 
la ^ b ^ c. Maintaining the prolate shape relies on dark 
matter staying on box orbits. A prime characteristics of box 
orbits is that they make close passages past the halo center. 
As baryons cool and collapse into the center, they deepen 
the potential well there and scatter material from box orbits 
into rounder loop or tube orbits. Tube orbits have an oblate 
shape, a ~ fo > c. Figs. 4 & 5 show the gradual transition of 
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Figure 4. Ratio of minor-to-major axes (c/a, left), middle-to-major axes {b/a, center) and triaxiality (T, right) as a fnnction of radius. DM-only simulations 
are depicted in black while hydro (NIHAO) simulations are depicted in red. Both the integral (solid line) and differential (dotted line) values are shown. 


inner halo shapes from prolate towards oblate as the mass 
increases: b/a increases to 1 more quickly than c/a. 

In order to better quantify the difference in shape be¬ 
tween DM-only and NIHAO simulations on a galaxy-by- 
galaxy base, we have decided to compare results for all three 


shape diagnostics {b/a, c/a and T) at three specific radii: 
i?vir, 12% i?vir and 5% i?vir®. 


® The difference in Avirbetween DM-only and Hydro is always of 
the order of few percent and hence does bias our results. 
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FigUrG 5. Same as Fig, 4 for all galaxies in the four mass bins M3-M0 (MO < Ml < M2 < M3 oj 10°°Mq). The shaded area represents the Icr scatter from 
galaxy to galaxy in the respective mass bin. Only the integral shape (solid line) is shown. 


The virial radius has been chosen to have the possibil¬ 
ity to make a direct comparison with previous studies based 
on DM-only simulations (Jing & Suto 2002; Allgood et al. 
2006; Bett et al. 2007; Maccio et al. 2008). On the other 
hand, observationally measuring the shape of dark matter 
haloes at the virial radius is a quite difficult task since there 


are very few (if any) tracers of the DM shape (or potential) 
that extend so far from the galaxy center. We have then de¬ 
cided to look at the halo shape at 0.12 x i?vir, which is in 
the middle of the 20-60 kpc radius range of the halo that 
Ibata et al. (2001) state shaped the Sagittarius stream. Fi¬ 
nally we also look at the shape in the very inner part of the 
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Figlir© 6. Axis ratios (c/a., b/a) and triaxiality parameter (T) for all the galaxies as a function of their virial (total) mass. The first row shows results at 
the virial radius, the middle row results at 0.12 of fiyir S'lid the last row at 5% of R-vir- R-ed diamonds represent the NIHAO galaxies, while black squares 
show the corresponding results for the DM-only runs. The solid lines in the first and second panel show the relation from Maccio et al. (2008). The green 
circle in the middle left panel shows the halo shape measurement for the Milky Way from Ibata et al. (2001). The red (black) line shows the median value 
for the hydro (DM) results. 


halo (5% -Rvir), which is the most prone to be affected by 
baryonic effects. 

Fig. 6 shows the value of the three shape diagnostics 
as a function of the total mass of the halo, at our three 
reference radii. At the virial radius, at all masses the NIHAO 
galaxies and their DM-only counterparts show very similar 
values for the axis ratio and Triaxiality, and in agreement 
with results obtained from larger samples of simulated dark 
matter haloes (Maccio et al. 2008, black solid line in the first 
and second panel). 

At smaller radii (second and third rows) there is a 
steady increase of c/a (and b/a) with Mvir in the NI¬ 
HAO simulations, which brings the simulated values in good 
agreement with the MW observations of Ibata et al. (2001, 
green circle). The last column in Fig. 6 shows the triaxiality 
parameter vs halo mass, confirming that CDM haloes are 
typically prolate (black squares). Interestingly, when galaxy 
formation is included the full range of halo triaxialities is 
possible (red diamonds). 

The difference between Hydro and DM-only simulations 
can be better appreciated in Fig. 7, where we show the ra¬ 
tio between the halo shapes c/a (left), b/a (center) and T 


(right) between collisional and collisionless simulations as a 
function of halo mass. The figure shows that there is a con¬ 
sistent shift for the inner halo shape from the DM-only to 
the baryon simulation and the trend of halo shapes becomes 
more spherical with increasing halo mass is clearly visible. 

For the left middle panel (0.12 f?vir, the most accessible 
to observations), we decided to fit such a ratio with a simple 
S-shape function: 

+ 1 

where M is the virial mass of the halo. We fixed the value 
of the S 2 parameter to 1.0 and fit for the other ones using 
the Levenberg-Marquardt method, results are reported in 
Table 2. The final fitting function is shown by the black line 
in Fig. 7, while the grey area represents the Icr = 0.156 
scatter around the mean. 

The trend of shape change with halo mass can be under¬ 
stood as a consequence of the increased efficiency of star for¬ 
mation of massive haloes (in our sample, e.g. 10^^ Mq ) with 
respect to low mass ones, as implied by abundance matching 
(Moster et al. 2010). Fig. 8 shows the different shape mea- 
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Figlir© T. Ratio between the halo shape (c/a) in the NIHAO and DM-only simulations as a function of the halo mass. The first row shows results at the 
virial radius, the middle row results at 0.12 of J?yij.and the last row at 5% of Rvir- The dashed line in the middle-left panel is the fitting function provided 
in Eq. 1, the grey area is the Icr scatter around the mean. The blue line shows the median value. 


Table 2. Fitting parameters describing the ratio between the halo 
shape, c/a, at 0.12 i?vir iri hydro and DM-only simulations 
(middle left panel of Fig. 7). 


Si S 2 Mo [Mq] /3 

1.848 1.0 3.1 X IQii 1.49 


surements versus the ratio between the stellar and the total 
mass. Since there are no stars in the DM-only simulation we 
use the empirical formula of Moster et al. (2010) to assign 
a stellar mass to each halo (black squares). 

As expected for a very low values of M*/Mvir < 0.01 the 
influence of baryons is minimal and the halo shape does not 
change substantially between DM-only and hydro simula¬ 
tions. For larger values of the stellar-to-total mass ratio, the 
two distributions tend to diverge, with the NIHAO galaxies 
showing a substantially rounder halo shape. 

A similar trend also applies to the Triaxiality param¬ 
eter (last columns). When galaxy formation is very ineffi¬ 
cient (Mstar/Mvir < 0.01) halo retain their prolate shape. 
Above (Mstar/A/vir > 0.01) haloes can become both triaxial 
T ~ 0.6 and close to oblate T ^ 0.2. This is quite impor¬ 


tant since, when embedded within a triaxial dark matter 
halo there can be systematic differences between the rota¬ 
tion speed and the circular velocity in gaseous discs, which 
has important implications for interpreting observations of 
dark matter density profiles (Hayashi & Navarro 2006). 

The shape of the halo leaves an imprint on the 2D kine¬ 
matics (Kuzio de Naray & Kaufmann 2011), with triaxial 
haloes yielding twists along the minor axis, and spherical 
(or aligned axisymmetric) haloes yielding symmetric veloc¬ 
ity fields. Our simulations thus predict a wider diversity of 
kinematic structures than would be inferred from DM-only 
simulations. 

Recently Kazantzidis et al. (2010) have used controlled 
(N-body) experiments to study the effect of the growth of 
a disc onto a triaxial halo. They found that the net effect 
depends weakly on the time scale of the disc assembly but 
strongly on the overall gravitational importance of the disc. 

In order to test their results on the importance of 
the baryonic mass (stellar and gaseous)^ contribution to 
the local potential in Fig. 9 we show the baryonic frac- 


^ Since we have fully cosmological simulations we have decided 
to use the global baryonic contribution and not just the disc one. 
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Figure 8. Same as Fig. 6 but this time as a function of the star formation efficiency of the halo, parameterized as Af*/Afvir- The values for Af*/iVfvir for 
the DM-only simulation are obtained using the relation from Moster et al. (2010). The red (black) line shows the median value for the hydro (DM) results. 


tion, [Mgas(-R) + Mt{R)\/Mtot{R), computed at R=5%-Rvir 
(green) and 12%_Rvir (blue), as a function of halo mass. Let 
us underline that our galaxies are very realistic stellar masses 
(see Fig. 1) and hence our baryonic fractions should be min¬ 
imally effected by a possible over-cooling. A comparison of 
this figure with Fig. 7 shows that a difference between DM- 
only and hydrodynamical simulations becomes clear only 
when the baryonic (mostly stellar) mass is at least of the 
order of 10% of the total one, supporting the findings of 
Kazantzidis et al. (2010). 

Finally before concluding our section on the DM halo 
shape we want to look at the effects of the shape of the 
stellar component (i.e. the galaxy morphology) on the final 
halo shape. Fig. 10 shows the relation between the shape 
of the stellar component (computed at 5% i?vir, radius that 
includes more than 90% of all stars) and the relative change 
in the inner shape (12% i?vir upper panel and 5% i?vir lower 
panel) between NIHAO and the DM-only run. In order to 
avoid to mix the effect of halo mass with the effect of stellar 
morphology this plot has been done using only galaxies in 
the M2 mass bin. 

On the larger spatial scale the correlation (if present) 
is very weak, while a positive trend is more apparent on the 
very small scale of 5% i?vir. On this scale it looks that a more 
spherical stellar distribution causes stronger change in the 
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Figure 9. Baryon fraction at 12% (green) and 5% (blue) of Hvir a 
function of total mass. It is clear that baryons start to be significant only 
above a total mass of few IO^^Mq. 


DM halo shape. This correlation can be explained by the fact 
that rounder stellar bodies (i.e. bulges) are more compact 
than disc ones and hence more effective in modifying the 
orbit of DM particles (e.g. Debattista et al. 2008). 
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lower panel) between NIHAO and the DM-only run. 
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3.2 Dark matter pseudo phase-space density 

Moving beyond the morphology of the dark matter, it is 
possible to also consider its kinematics. Taylor & Navarro 
(2001) defined “pseudo phase space density” as a simple re¬ 
lationship between matter density and velocity dispersion, 
a quantity that describes the matter distribution and kine¬ 
matics of the dark matter together. They defined pseudo 
phase space density as: 

Q{r) = p{r)/a^{r), 

where a{r),p{r) are the velocity dispersion and density of 
the halo, respectively. Taylor & Navarro (2001) found that 
this simple combination of properties serves as a useful probe 
for understanding the origin of the universal DM halo pro¬ 
files. 

Using DM-only simulations (Taylor & Navarro 2001) 
found that the pseudo phase space density follows a sim¬ 
ple power law, Q{r) oc r^, with y ~ —1.875. Since 
there is mounting evidence that baryons modify the 
DM density profile, p{r), (e.g., Mashchenko et al. 2006; 
Governato et al. 2010; Teyssier et al. 2012; Di Cintio et al. 
2014a; Toilet et al. 2016), it is worth checking whether 
baryons also reshape the Q profile. 

Fig. 11 presents a comparison between the matter den¬ 
sity profile (right) and the pseudo phase-space density Q 
profile (left) for our “test” galaxies (the same galaxies shown 
in Fig. 4). Each density profile includes three lines: the den¬ 
sity in the DM-only simulations (black solid line), the DM 
density in the hydro simulations (red solid line) and the 
total density profile (dark-|-stars-|-gas) in the hydro simula¬ 
tions (red dotted line). While profiles from DM-only simula¬ 
tions have universal Einasto-like shapes (e.g. Navarro et al. 
2004; Merritt et al. 2005; Dutton & Maccio 2014), profiles 
from hydro simulations exhibit a core for low mass haloes 
which gradually steepens with halo mass, becoming even 
more cuspy DM-only simulations of the most massive galax¬ 
ies in our sample (see Toilet et al. 2016 for a thorough dis¬ 
cussion of the modification of the DM density profiles in the 
NIHAO simulations). 

The Q profiles show a different behavior. Q in hydro 
simulations is always lower than its N-body counterpart. 
The profiles also flatten in the center even when the density 


profiles do not. The flattening of Q in the lowest mass halo 
(bottom panel) reflects the flattening of the matter density, 
which decreases the numerator in the definition of Q. For 
the other three panels, there is no such discrepancy in the 
matter density profile. However, the “total” density profile 
in the hydro simulation does get steeper in all three cases. 
The deeper global potential well causes the velocity disper¬ 
sion of the DM a to increase towards the center, thus flat¬ 
tening the inner part of the Q radial profile. As a conse¬ 
quence the pseudo phase-space density profile of the Dark 
Matter component departs from simple power law behavior 
at practically all mass scales probed by our hydrodynamical 
simulations. 

A comparison of Fig. 4 and Fig. 11 also suggests that 
the violent gas motions caused by baryons in the center 
of low mass galaxies which are believed to be responsi¬ 
ble of the flattening of the central dark matter cusps (e.g. 
Pontzen & Governato (2012)), also works to make the dark 
matter rounder. In the low mass g4.99el0, baryons make 
little change to the outer shape of the dark matter. How¬ 
ever, the dark matter is almost completely spherical inside 
2 kpc, the same region in which the density profile has been 
flattened. 

Fig. 12 shows the same quantities as in Fig. 11 but 
this time averaged for all galaxies in the different mass bins 
M0-M3. The averaged profiles confirm the trends already 
seen for the four “test” galaxies. In the lowest mass bin 
(MO) there is quite small difference between the averaged Q 
computed in the DM-only or NIHAO simulations, suggesting 
that stochasticity of star formation is quite strong on those 
mass scales as already pointed out in works dealing with 
the flattening of the density profiles (e.g. Onorbe et al. 2015; 
Toilet et al. 2016) 

In Fig. 13 we try to summarize our findings by showing 
the ratio of the values of Q at two different radii in NIHAO 
and in DM-only simulations: 0.05 i?virand 0.12 Rvii- At low 
masses (M« 10^° Mq ) there is no much difference between 
the two Qs, at higher masses the hydro simulations show 
on average a lower value of Q, in agreement with findings 
based on single profiles. Most likely this this strong depen¬ 
dence of the Q ratio on halo mass is due to the term in 
the denominator of the definition of the pseudo phase-space 
density. The density when computed at 5% or i?vir (lower 
panel of figure 13) shows an average lower value in the 
hydro simulations up to a mass of few 10^^ solar masses, 
and similar values above this mass, consistent with results 
on halo expansion and contraction reported in Toilet et al. 
(2016). 

3.3 Dark matter velocity distribution 

As outlined in the introduction, the velocity distribution is 
crucial for the detection of DM in the Milky-Way. For this 
reason together with the results of our test galaxies, and 
the corresponding averaged mass bins (M0-M3), we will also 
show results for four more single galaxies namely: g8.26ell, 
gl.l2el2, gl.77el2, g2.79el2. 

These galaxies have been selected in order to have a 
stellar and halo mass similar to the one of our own Galaxy 
(« IO^^Mq), three of them (g8.26ell, 1.77el2 and g2.79el2) 
are strongly disc dominated with a disc to total ration of 
D/T > 0.6 (computed according to Obreja et al. (2016)), 
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Figlir© 11. The density (left) and pseudo phase space density (right) of DM-only (black) and NIHAO (red) galaxies for our four test objects. Solid lines 
represent the density and Q values of the dark matter particles, while the red dotted line represents the total density of the star, gas, and dark matter 
particles. 
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Figure 12. The density (left) and pseudo phase space density (right) of DM-only (black) and NIHAO (red) galaxies averaged in the for mass bins M0-M3 
(MO < Ml < M2 < M3 ~ 10^^ Mq). The grey area represents the 1 <t scatter around the mean. 
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Figure 13. The top panel shows the ratio of the values of Q in the in 
dissipationless and hydro simulations computed at 0.05 i?vir(0-12 
the middle panel). The outliers in the plot are due to merging galaxies 
(haloes) which are at a sligtly different merging state in the Hydro and DM 
simulations. The bottom panel shows the values of p at 5% of fhe 

DM (black) and the NIHAO (red) simulations. 


while gl.l2el2 is quite compact with the disc accounting for 
~ 10% of the total stellar mass. Edge-on and face-on images 
of these four galaxies are shown in figure 15. Each image 
is 50 kpc on a side and was created using the Monte Carlo 
radiative transfer code sunrise (Jonsson 2006). 

Fig. 14 shows the velocity distribution of dark matter 
on two different scales for our test galaxies. The first row 
of Fig. 14 shows the distribution of all DM particles within 
the virial radius while the second row shows the distribu¬ 
tion only inside the solar neighborhood. This quantity is 


defined as a shell of radius between 7 kpc < r < 9 kpc for 
galaxies in the M3 mass bin, and it is then rescaled to a sim¬ 
ilar fraction of the virial radius (around 4%) for lower mass 
galaxies. In principle the solar neighborhood should be de¬ 
fined as a ring in the plane of the disc with radius between 
7 and 9 kpc. We tested on our more disc dominated galaxies 
(g8.26ell, gl.77el2, gl.92el2, g2.79el2) that the results do 
not change at all if we use a spherical shell instead, which 
has two advantages: it increases the number of DM particles 
(and thus reduces the numerical noise) and can be applied 
to both hydro and DM-only simulations. The similarity of 
results for a shell and a ring suggests that the DM distribu¬ 
tion is spherically symmetric and no “dark” disc is present 
in our simulations. The black lines in figure 14 show results 
from DM-only simulations while red lines are for the NIHAO 
galaxies. At the virial radius, in both types of simulations, a 
Maxwellian distribution well represents the global velocity 
distribution: 


f{x) = mi 


^2g-xV(2m|) 

ml 


( 2 ) 


The best fit Maxwell functions are shown as dashed lines in 
Fig. 14. 

When restricted to the solar neighborhood, the DM- 
only and the hydro velocity distributions differ substantially 
in all galaxies. The DM-only simulations can still be well fit 
with a Maxwellian distribution, (even if they show a slightly 
larger tail at high velocity). Kuhlen et al. (2010) found sim¬ 
ilar agreement using higher resolution simulations. In the 
hydro case, the velocity distribution for the most massive 
test galaxy (gl.92el2) is much more symmetric around the 
maximum value, which increases substantially, and then the 
distribution falls quite rapidly at high velocities. This differ¬ 
ence is most likely due to the quite significant halo contrac¬ 
tion for this galaxy (see Fig. 11) which boosts the local DM 
velocity dispersion. 

This strong variation at the very high mass end of the 
NIHAO galaxies is also confirmed by a closed inspection of 
the velocity distribution of the four additional high mass 
galaxies shown in Fig. 16. With increasing mass the distri¬ 
bution becomes more Gaussian and the peak moves towards 
larger velocities, but the fall-off at the high velocity end 
is quicker than in the DM-only case. Despite the different 
aspect (see figure 15), the four galaxies have quite similar 
behavior, suggesting that morphology is not the main driver 
of the changes in the velocity distribution, which seems to 
be more related to the total stellar (or halo) mass. 

For smaller masses (g5.02ell, gl.OSell, g4.99el0) the 
peak of the distribution moves towards lower velocities, 
possibly related to the different halo response on these 
mass scales, where either expansion or no reaction (for 
M*/Mhaio < 10“'") is expected (Toilet et al. 2016). 

As a consequence of baryonic effects, the local velocity 
distribution of MW-like galaxies is better fit by a Gaussian 
distribution: 

= (3) 


The Gaussian fit is shown in the lower panels of Figs. 14 & 16 
by the (red) dotted line that clearly provides a better fit than 
a Maxwellian (red) dashed line. The fitting parameters for 
both distributions (Gaussian and Maxwellian) are reported 
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Figure 14. The dark matter particles velocity distribution in our four test galaxies. The dashed lines indicate a Maxwellian fit to the simulated curves. 
The top panels show the global velocity distribution, while the bottom panels show local measurements taken at the solar position ( 7 kpc < r < 9 kpc) for 
galaxies in the M3 mass bin, and a similar fraction (ss 5%) of the virial radius for lower mass galaxies. For the local velocity distribution the dotted line 
shows a Gaussian fit to the NIHAO curve. 



Figlir© 15. Face-on (upper panels) and edge-on (lower panels) views of galaxies: g8.26ell, gl.l2el2, gl.77el2, g2.79el2 after processing through the Monte 
Carlo radiative transfer code SUNRISE. Images are 50 kpc on a side. 


in Table 3, all fits have been performed using the Levenberg 
and Marquardt algorithm 

Finally Fig. 17 shows the average velocity distribution 
and the virial radius (upper panels) and at the solar radius 
(lower panels). While the trend seen is individual galaxies 
confirmed in the M0-M2 bins (with gl.99el0 being a bit bor¬ 
derline), this figure also illustrates the non negligible scatter 
in the galaxy-to-galaxy variation for the local velocity dis¬ 
tribution. 

In the highest mass bin (M3) the results previously 
highlighted get washed away due to the large scatter galaxy- 
by-galaxy as marked by the extended red area. This scatter 
is due to the rapid change in halo response (from expan¬ 


sion to contraction) across this mass bin. For lower masses 
the average values confirm trends seen previously on single 
galaxies, namely more symmetric distributions and a lower 
tail at high velocities. 

On the Milky-Way scale our findings confirm ear¬ 
lier results based on the single halo ERIS simulation 
(Guedes et al. 2011) and presented by Pillepich et al. 
(2014). In their paper Pillepich and collaborators also re¬ 
ported a suppression of the wings and a more symmetric 
shape for the velocity distribution. In our simulations the 
differences between the DM-only case and the hydro simula¬ 
tions seem to be even larger. This discrepancy might be due 
to the stronger feedback model (which leads to a stronger 
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Figure 16. The dark matter particles velocity distribution in four additional galaxies with stellar and DM masses similar to the Milky Way: g8.26ell, 
gl.l2el2, gl.77el2, g2.79el2. Symbols and lines are the same as in Fig. 14. 
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Figur© IT. Same as Fig. 14 but averaged in the four mass bins M0-M3 (MO < Ml < M2 < M3 ~ IO^^Mq). 


baryonic impact) we have adopted in our simulations in or¬ 
der to balance the metallicity dependent gas cooling, which 
was ignored in the original ERIS simulation. 

As extensively discussed in Pillepich et ah, the suppres¬ 
sion of the tail of the distribution at high velocities has im¬ 
portant consequences on the interpretation and the compar¬ 
ison of different direct detection experiments. For example 
it relaxes the tension between the (possible) signal of dark 
matter scattering reported by CDMS-Si (Agnese & et al. 
2013) and the exclusion of such a signal from the Xenon- 
100 experiment (Aprile & et al. 2012). We refer the reader 
to Mao et al. (2014) and Pillepich et al. (2014) for a more 
thorough discussion. 


4 CONCLUSIONS 

We used the NIHAO simulation suite (Wang et al. 2015) to 
investigate the impact of galaxy formation on the properties 
of the dark matter distribution within haloes. 

The NIHAO suite is a large simulation campaign aim¬ 
ing to produce a large sample of high resolution simulated 
galaxies in a cosmological context. It is an extension of the 


MaGICC simulations (Stinson et al. 2013) and it has been 
performed with an improved version of the SPH GASOLINE 
code (Keller et al. 2014) which fixes the well known prob¬ 
lems of particle based hydrodynamical codes (Agertz et al. 
2007). The NIHAO project counts more than 90 simu¬ 
lated galaxies across two orders of magnitude in halo mass 
(10^'’ — 10^^ Mq), with each of the galaxies resolved with at 
least 4 X 10® elements. 

The NIHAO galaxies have been very successful in repro¬ 
ducing stellar to halo mass ratio on more than five orders of 
magnitude in stellar mass: from 10® to IO^^Mq. They also 
show very realistic star formation histories for their stellar 
and halo masses. Finally NIHAO galaxies are also consistent 
with the observed gas content of galaxy discs and haloes 
(Stinson et al. 2015; Wang et al. 2016; Gutcke et al. 2016). 
Thanks to the unprecedented combination of high resolution 
and large statistical sample, the NIHAO suite offers a unique 
set of objects to study the distribution response of several 
DM properties to galaxy formation. Moreover since for each 
galaxy we have, at the same resolution, an N-body only (col¬ 
lisionless) simulation and a full hydrodynamical simulation. 
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Table 3. Parameters describing the velocity distribution function 
at the virial radius (first two rows) and at the solar radius (second 
two rows) for the five Milky-Way like galaxies. 


g8.26ell 

mi 

m2 

91 

9 

cr 

Rvir (dm) 

0.789 

109.5 

0.0055 

158.9 

70.42 

Avir (hydro) 

0.794 

105.7 

0.0053 

164.3 

70.86 

Rq (dm) 

0.776 

131.0 

0.0044 

195.2 

103.37 

Rq (hydro) 

0.903 

135.5 

0.0061 

211.8 

64.69 

gl.l2el2 

mi 

m2 

91 

9 

cr 

(dm) 

0.787 

101.9 

0.0055 

155.6 

72.28 

Avir (hydro) 

0.755 

110.7 

0.0049 

165.0 

76.96 

Rq (dm) 

0.754 

127.0 

0.0041 

187.3 

104.2 

Rq (hydro) 

0.905 

147.7 

0.0066 

241.9 

46.79 

gl.77el2 

mi 

m2 

91 

9 

cr 

i?vir (dm) 

0.725 

129.8 

0.0042 

195.5 

99.00 

Avir (hydro) 

0.713 

137.0 

0.0039 

206.6 

99.80 

Rq (dm) 

0.712 

137.1 

0.0032 

189.2 

115.22 

Rq (hydro) 

0.929 

208.3 

0.0047 

328.7 

76.69 

gl.92el2 

mi 

m2 

91 

9 

cr 

i?vir (dm) 

0.739 

144.5 

0.0035 

211.6 

122.3 

Avir (hydro) 

0.699 

140.37 

0.0033 

203.6 

119.5 

Rq (dm) 

0.735 

159.5 

0.0030 

232.0 

115.2 

Rq (hydro) 

0.969 

218.0 

0.0058 

348.1 

60.49 

g2.79el2 

mi 

m2 

91 

9 

cr 

i?vir (dm) 

0.783 

155.0 

0.0036 

231.9 

110.0 

Avir (hydro) 

0.774 

160.1 

0.0035 

239.4 

109.2 

Rq (dm) 

0.722 

170.9 

0.0027 

249.0 

133.8 

Rq (hydro) 

0.870 

256.6 

0.0045 

401.5 

81.91 


we are able to assess the effect of baryons on a halo-by-halo 
basis. 

In this study , we focused on three key properties of 
dark matter haloes: the halo shape, the radial profile of the 
pseudo phase-space density and the dark matter velocity 
distribution both global and in the solar neighborhood. Our 
results can be summarized as follows: 

(i) The shape of the dark matter halo within the virial 
radius is similar between DM-only and hydro simulations. 
At smaller radii, however, the hydro simulations become 
rounder. There is a strong mass dependence to the differ¬ 
ence between the inner halo shape (measured at 12 % of the 
virial radius) from the DM-only simulations and hydro sim¬ 
ulations. At low masses (< lO^^M©) the dark matter halo 
tends to retain its original triaxial shape, while at higher 
masses (« 1 O^^M 0 ) the inner halo becomes more spherical 
with an average minor to major axis ratio (c/a) of 0.8 . This 
brings numerical predictions into good agreement with esti¬ 
mates of the inner halo shape in our own Galaxy. We show 
that the mass dependence of the variation of the halo shape 
is related to the increase of star formation efficiency with 
halo mass, which raises the contribution of stars and gas to 
the overall potential. 

We provide a simple fitting formula that relates the 
change in the axis ratio c/a with the ratio between stellar 
mass and halo mass, and in principle allows the prediction 


of the shape of a galaxy DM halo from its stellar and 
dynamical mass. 

(ii) In hydrodynamical simulations the radial behavior 
of the dark matter pseudo phase space-density Q = p/cr^ is 
not always well represented by a single power law. At total 
masses M /i 10^^ Mq the Q radial profile shows a flattening 
towards the center of the halo. This is related to the change 
in the DM density profile which strongly departs from the 
pure DM results in the NIHAO simulations (Toilet et al. 
2016). Overall hydro simulations have a lower value of Q 
compared to the Nbody case, when it is measured at a fix 
fraction of the virial radius. 

(iii) The velocity distribution of the dark matter particles 
within the virial radius in the hydro simulations is still well 
represented by a Maxwellian distribution, and it is similar 
to the DM-only case at all mass scales. 

When we restrict our analysis to the the solar neighbor¬ 
hood (7 kpc < r < 9 kpc) we find that in the hydro sim¬ 
ulations the velocity distribution functional form strongly 
depends on the halo mass. 

At low halo masses {M ~ 10^° Mq the Hydro and DM- 
only simulations show a similar behavior, when we move 
to higher masses {M ~ lO^^M©) the velocity distribution 
becomes progressively more symmetric and the velocty peak 
moves towards lower values w.r.t the Nbody case, then in our 
most massive bin (M ~ 10 ^^Mq) the distribution is again 
in agreement with the Nbody case. We tentatively ascribed 
this trend to the different reaction of the DM distribution 
as function of increasing halo mass from few 10 ® Mq to 
M> 10^^ Mq (namely no effect, halo expansion, no effect 
and halo contraction) as described in details in Toilet et al. 
(2016) and Dutton et al. (2016a). 

To better study this effect, we isolated five galaxies that 
for stellar and total mass resemble our own Milky Way. For 
these galaxies the maximum of the velocity distribution 
in the hydro simulations moves to higher velocities w.r.t 
the Nbody case, due to an overall halo contraction in 
these galaxies, and we find very little correlation with the 
galaxy morphology. In our Milky Way analogs the velocity 
distribution is well fitted by a Gaussian and we provide 
the fitting parameters of the distribution for five different 
galaxies. We also stress that the lack of high velocity parti¬ 
cles has important consequences for the interpretation and 
comparison of Dark Matter direct detection experiments. 

Our results show that baryons have important effects on 
the dark matter not only in the very inner part (e.g. leading 
to expansion and contraction) but also on the global proper¬ 
ties of dark matter. The understanding of the nature of dark 
matter and the comparison of theoretical predictions with 
observational data can no longer rely on pure collisionless 
simulations, but must include the effects of visible matter. 
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